/*
 * EBFV1_loadPreprocessor.cpp
 *
 *  Created on: 16/09/2012
 *      Author: rogsoares
 */

#include "EBFV1_preprocessor.h"

namespace PRS{
	void createVertexMapping(TMesh *pTMesh_data);
	void definesEdgesVersorsAndLengths(TMesh* pTMesh_data);
	void getAllVertices(TMesh* pTMesh_data);
	void getAllEdges(TMesh* pTMesh_data);

	void EBFV1_loadPreProcessorData(string fname, TMesh *pTMesh_data){
		if (!0){
			cout << "loading pre-processor data begin.\n";
		}
		char filename[256];
		sprintf(filename,"%s-%d-of-%d.dat",fname.c_str(),0,1);
		PetscSynchronizedPrintf(PETSC_COMM_WORLD,"pre-processor file: %s\n",filename);
		PetscSynchronizedFlush(PETSC_COMM_WORLD);
		MPI_Barrier(PETSC_COMM_WORLD);

		ifstream fid;
		fid.open(filename);

		if (!fid.is_open()){
			char str[256];
			sprintf(str,"File %s could not be opened or it does not exist\n",filename);
			throw Exception(__LINE__,__FILE__,str);
		}

		loadNodalVolumes(fid,pTMesh_data);
		//printVerticesData(pTMesh_data);STOP();

		loadCij(fid,pTMesh_data);
		//printEdgesData(pTMesh_data);

		loadDij(fid,pTMesh_data);
		//printBdryElemData(pTMesh_data);

		// Number of global edges:
		string line;
		getline(fid,line,'\n');
		getline(fid,line,'\n');
		getline(fid,line,'\n');
		pTMesh_data->setNumGEdges( atoi(line.c_str()) );
		getline(fid,line,'\n');
		getline(fid,line,'\n');
		pTMesh_data->setNumLocalEdges(atoi(line.c_str()) );
		fid.close();

		/*
		 * Edges need vertices data. A vertex ID mapping is necessary.
		 */
		createVertexMapping(pTMesh_data);
		definesEdgesVersorsAndLengths(pTMesh_data);
		getAllVertices(pTMesh_data);
		getAllEdges(pTMesh_data);
	}

	void createVertexMapping(TMesh *pTMesh_data){
		cout << "Starting mapping... ";
		std::map<int,VertexData*> map;
		VertexData *vdata;
		int nrows, ncols;

		for(int i=0; i<pTMesh_data->getNumDomains(); i++){
			pTMesh_data->vertices_list[i].getSize(nrows,ncols);
			for(int j=0; j<nrows; j++){
				vdata = pTMesh_data->vertices_list[i].getValue(j,0);
				map[vdata->ID] = vdata;
			}

			EdgeData *edata;
			pTMesh_data->edge_list[i].getSize(nrows,ncols);
			for(int j=0; j<nrows; j++){
				edata = pTMesh_data->edge_list[i].getValue(j,0);
				edata->vertex_0 = map[edata->ID_0];
				edata->vertex_1 = map[edata->ID_1];
			}

			BdryElementData *bedata;
			int dim = pTMesh_data->getMeshDim();
			pTMesh_data->bdry_elem_list[i].getSize(nrows,ncols);
			for(int j=0; j<nrows; j++){
				bedata = pTMesh_data->bdry_elem_list[i].getValue(j,0);
				bedata->vertex_0 = map[bedata->ID_0];
				bedata->vertex_1 = map[bedata->ID_1];
				if(dim==3){
					bedata->vertex_2 = map[bedata->ID_2];
				}
			}
			map.clear();
		}
		cout << "done!";

		// verify what has been done
	//	for(int i=0; i<pTMesh_data->getNumDomains(); i++){
	//				EdgeData *edata;
	//				pTMesh_data->edge_list[i].getSize(nrows,ncols);
	//				for(int j=0; j<nrows; j++){
	//					edata = pTMesh_data->edge_list[i].getValue(j,0);
	//					cout << edata->vertex_0->ID << "\t" << edata->vertex_0->volume << "\t" << edata->vertex_0->numRC << "\t";
	//					cout << edata->vertex_1->ID << "\t" << edata->vertex_1->volume << "\t" << edata->vertex_1->numRC << "\n";
	//				}
	//
	//		BdryElementData *bedata;
	//		int dim = pTMesh_data->getMeshDim();
	//		pTMesh_data->bdry_elem_list[i].getSize(nrows,ncols);
	//		for(int j=0; j<nrows; j++){
	//			bedata = pTMesh_data->bdry_elem_list[i].getValue(j,0);
	//			cout << bedata->vertex_0->ID << "\t" << bedata->vertex_0->volume << "\t" << bedata->vertex_0->numRC << "\t";
	//			cout << bedata->vertex_1->ID << "\t" << bedata->vertex_1->volume << "\t" << bedata->vertex_1->numRC << "\n";
	//		}
	//	}
	}

	void definesEdgesVersorsAndLengths(TMesh* pTMesh_data){
		double minEdgeLength = 1e10;
		int nrows, ncols;
		int dim = pTMesh_data->getMeshDim();
		int ndom = pTMesh_data->getNumDomains();
		for(int i=0; i<ndom; i++){

			// Domain's edges
			pTMesh_data->edge_list[i].getSize(nrows,ncols);
			for(int j=0; j<nrows; j++){
				EdgeData *edata = pTMesh_data->edge_list[i].getValue(j,0);
				VertexData* vdata[2] = {edata->vertex_0, edata->vertex_1};
				double IJ_vector[3] = {vdata[1]->x - vdata[0]->x,vdata[1]->y - vdata[0]->y, vdata[1]->z - vdata[0]->z};
				double sign = (edata->vertex_0->ID > edata->vertex_1->ID)?-1.:1.;
				edata->length = norm_L2(IJ_vector,dim);
				edata->versor = new double[dim];
				for(int k=0; k<dim; k++){
					edata->versor[k] = sign*IJ_vector[k]/edata->length;
				}
				minEdgeLength = std::min(minEdgeLength,edata->length);
				//cout << "edata->length: " << edata->length << "\tmax length: " << minEdgeLength << endl;
			}

			if (dim==2){
				// Boundaries' edges
				pTMesh_data->bdry_elem_list[i].getSize(nrows,ncols);
				for(int j=0; j<nrows; j++){
					BdryElementData *edata = pTMesh_data->bdry_elem_list[i].getValue(j,0);
					VertexData* vdata[2] = {edata->vertex_0, edata->vertex_1};
					double IJ_vector[3] = {vdata[1]->x - vdata[0]->x,vdata[1]->y - vdata[0]->y, vdata[1]->z - vdata[0]->z};
					double sign = (edata->vertex_0->ID > edata->vertex_1->ID)?-1.:1.;
					double length = norm_L2(IJ_vector,dim);
					edata->versor = new double[dim];
					for(int k=0; k<dim; k++){
						edata->versor[k] = sign*IJ_vector[k]/length;
					}
				}
			}
		}


		// todo: broadcast maximum edge length for all processes
		pTMesh_data->setMaximumEdgeLentgh(minEdgeLength);
	}

	void getAllVertices(TMesh* pTMesh_data){
		int nrows, ncols, ndom;
		std::map<int,VertexData*> mapAllVertices;
		ndom = pTMesh_data->getNumDomains();
		for(int i=0; i<ndom; i++){
			pTMesh_data->vertices_list[i].getSize(nrows,ncols);
			for(int j=0; j<nrows; j++){
				VertexData* vdata = pTMesh_data->vertices_list[i].getValue(j,0);
				mapAllVertices[vdata->ID] = vdata;
			}
		}
		int size = (int)mapAllVertices.size();
		int row = 0;
		int ID_local = 0;	// for visualization purposes only
		pTMesh_data->allVertices.allocateMemory(size,1);
		for (std::map<int,VertexData*>::iterator iter = mapAllVertices.begin(); iter!=mapAllVertices.end(); iter++){
			pTMesh_data->allVertices.setValue(row,0,iter->second);
			VertexData* vdata = pTMesh_data->allVertices.getValue(row,0);
			vdata->ID_local = ID_local;
			ID_local++;
			row++;
		}
		mapAllVertices.clear();
	}

	// for visualization purposes only
	// todo: not cjeked for multidomains
	void getAllEdges(TMesh* pTMesh_data){
		int nrows, ncols, ndom, i, j;
		std::set<int> bdryID;
		std::set<BdryElementData*> bdryEdges;

		// get ID from all nodes laying on internal boundaries (between sub-domains) only
		ndom = pTMesh_data->getNumDomains();
		for(i=0; i<ndom; i++){
			pTMesh_data->bdry_elem_list[i].getSize(nrows,ncols);
			for(j=0; j<nrows; j++){
				BdryElementData* bedata = pTMesh_data->bdry_elem_list[i].getValue(j,0);
				// look for edges between sub-domains
				if (bedata->dom1 && bedata->dom2){
					bdryID.insert( bedata->vertex_0->ID );
					bdryID.insert( bedata->vertex_1->ID );
					bdryEdges.insert(bedata);
				}
			}
		}

		// put all edges from all sub-domains into the same bag (edges between sub-domains will be repeated)
		std::set<EdgeData*> edgeBag;
		std::set<int>::iterator iter;
		for(i=0; i<ndom; i++){
			pTMesh_data->edge_list[i].getSize(nrows,ncols);
			for(j=0; j<nrows; j++){
				EdgeData* edata = pTMesh_data->edge_list[i].getValue(j,0);
				// insert only edges whose vertices IDs are not present into the bdryID container
				iter = bdryID.find(edata->vertex_0->ID);
				if (iter==bdryID.end()){
					iter = bdryID.find(edata->vertex_1->ID);
					if (iter==bdryID.end()){
						edgeBag.insert(edata);
					}
				}
			}
		}
//		cout << "nrows = " << nrows << endl;
//		cout << "edgeBag.size = " << edgeBag.size() << endl;STOP();

		// all edges matrix
		nrows = (int)edgeBag.size();
		nrows += (int)bdryEdges.size();
		cout << "nrows = " << nrows << endl;
		pTMesh_data->allEdges.allocateMemory(nrows,2);
		j = 0;
		for(std::set<EdgeData*>::iterator edataIter = edgeBag.begin(); edataIter!=edgeBag.end(); edataIter++){
			EdgeData* edata = *edataIter;
			pTMesh_data->allEdges.setValue(j,0,edata->vertex_0);
			pTMesh_data->allEdges.setValue(j,1,edata->vertex_1);
			j++;
		}
		for(std::set<BdryElementData*>::iterator bedataIter = bdryEdges.begin(); bedataIter!=bdryEdges.end(); bedataIter++){
			BdryElementData* edata = *bedataIter;
			pTMesh_data->allEdges.setValue(j,0,edata->vertex_0);
			pTMesh_data->allEdges.setValue(j,1,edata->vertex_1);
			j++;
		}

		bdryEdges.clear();
		edgeBag.clear();
		bdryID.clear();

		// let's subtract thoses edges which are between two sub-domains
		// A loop over all boundary edges will be necessary, but only those between two sub-domains will be used
		// External boundary edges are skipped
		std::set<EdgeData*> repeatedEdges;
	}
}
